clear
rng default
load IdSetM_final_x1
IdSetM{1}=IdSetM_final_x1;
load IdSetM_final_x2
IdSetM{2}=IdSetM_final_x2;
load IdSetM_final_x3
IdSetM{3}=IdSetM_final_x3; 
load IdSetM_final_x4
IdSetM{4}=IdSetM_final_x4;
load IdSetM_final_x5
IdSetM{5}=IdSetM_final_x5;
load IdSetW_final_y1
IdSetW{1}=IdSetW_final_y1;
load IdSetW_final_y2
IdSetW{2}=IdSetW_final_y2;
load IdSetW_final_y3
IdSetW{3}=IdSetW_final_y3; 
load IdSetW_final_y4
IdSetW{4}=IdSetW_final_y4;
load IdSetW_final_y5
IdSetW{5}=IdSetW_final_y5;

load U_CS
load V_CS
load PYX_cond
load PXY_cond

ntypes=4;

C_interval_m=cell(ntypes,ntypes);
C_interval_w=cell(ntypes,ntypes);
for i=1:ntypes-1
    for k=i+1:ntypes
        %Men
        U_i=IdSetM{i}(:,2:end);
        U_k=IdSetM{k}(:,2:end);
        C_tempk=sum(PYX_cond(k,1:end-1).*U_k,2);
        C_tempi=sum(PYX_cond(i,1:end-1).*U_i,2);
        C_interval_m{k,i}=[min(C_tempk)-max(C_tempi) max(C_tempk)-min(C_tempi)];
        
        %Women
        V_i=IdSetW{i}(:,2:end);
        V_k=IdSetW{k}(:,2:end);
        C_tempk=sum(PXY_cond(k,1:end-1).*V_k,2);
        C_tempi=sum(PXY_cond(i,1:end-1).*V_i,2);
        C_interval_w{k,i}=[min(C_tempk)-max(C_tempi) max(C_tempk)-min(C_tempi)];
    end
end


C_CS_m=zeros(ntypes,ntypes);
C_CS_w=zeros(ntypes,ntypes);
for i=1:ntypes-1
    for k=i+1:ntypes
        %Men
        U_i=U_CS(i,:);
        U_k=U_CS(k,:);
        C_CS_m(k,i)=sum(PYX_cond(k,1:end-1).*U_k,2)-sum(PYX_cond(i,1:end-1).*U_i,2);
        
        
        %Women
        V_i=V_CS(:,i).';
        V_k=V_CS(:,k).';
        C_CS_w(k,i)=sum(PXY_cond(k,1:end-1).*V_k,2)-sum(PXY_cond(i,1:end-1).*V_i,2);
    end
end


save('C_interval.mat', 'C_interval_m', 'C_interval_w')
save('C_CS.mat', 'C_CS_m', 'C_CS_w')
